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Summary 


A direct procedure is presented for locally bicubic interpolation on a structured, curvi- 
linear, two-dimensional grid. The physical (Cartesian) space is transformed to a compu- 
tational space in which the grid is uniform and rectangular by a generalized curvilinear 
coordinate transformation. Required partial derivative information is obtained by finite 
differences in the computational space. The partial derivatives in physical space are de- 
termined by repeated application of the chain rule for partial differentiation. A bilinear 
transformation is used to analytically transform the individual quadrilateral cells in phys- 
ical space into unit squares. The interpolation is performed within each unit square using 
a piecewise bicubic spline. 


1. Introduction 

Problems in computational fluid dynamics (CFD) are often solved using a numeri- 
cally generated grid which is structured but nonuniform, consisting of quadrilaterals in 
two dimensions and hexahedra in three dimensions. The grid is structured in that each 
cell contains an implied set of coordinate directions. A uniform rectangular grid in com- 
putational space is obtained from the physical space grid by a generalized curvilinear 
coordinate transformation (ref. 1). This facilitates the application of finite differences and 
finite volume methods. 

In many applications of CFD, such as multiblock algorithms using overlapping grids, 
surface definition, and graphical representation of flow fields, data which are known at grid 
points must be interpolated at intermediate points. In the overlapping grid technique de- 
scribed by Benek et al. (ref. 2), interpolation is accomplished using an iterative technique. 
The purpose of the present paper is to present an alternative approach which involves no 
iteration. 

The interpolation procedure presented uses a piecewise bicubic interpolant rather 
than the simpler bilinear interpolant. The bilinear interpolant has the often desirable 
property that it introduces no new extrema. However, the bicubic interpolant is more 
accurate and produces a smooth interpolant. Furthermore, bicubic interpolation is capable 
of preserving derivative information generated by a finite-difference flow solver. Important 
considerations regarding monotonicity and conservative properties of the interpolant are 
not addressed in this paper. Examples of alternative techniques for bivariate interpolation 
on a rectangular grid are given by Carlson and Fritsch (ref. 3) and Schiess (ref. 4). 

The generalized curvilinear coordinate transformation is given in Section 2, and the 
iterative interpolation technique of Benek et al. is reviewed in Section 3. The direct 
interpolation procedure is presented in Section 4. The Appendix contains the formulas for 
bicubic spline interpolation on a unit square. 
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2. Generalized Curvilinear Coordinate Transformation 

The generalized curvilinear coordinate transformation maps a nonuniform, nonrect- 
angular grid in physical (Cartesian) space to a uniform, rectangular computational grid. 
The Cartesian coordinates (x,t/) are mapped to the curvilinear coordinates (£, 77 ) by the 
following general transformation: 


v = v 

where we consider two space dimensions only and ignore time dependence for simplicity. 
The mapping is chosen such that the resulting grid in computational space is uniform, 
rectangular, and of unit length, i.e., 


= Arj = 1. (2) 

Hence, finite-difference representations of and d v are easily formulated. The mapping 
is one-to-one except at topological singularities or cuts, where one physical point may be 
mapped to many computational points. 

The chain rule for partial differentiation gives, in matrix form, 


d x 


'U 

Vx 


ki 

kJ 


Uv 

V y. 


kJ 


The coordinate transformation is generally not known analytically and hence must be 
determined numerically. When the roles of the coordinate systems in equation (3) are 
reversed, the chain rule gives 


k" 


*€ 

Vi 

d x 

kJ 
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Therefore, we must have 
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Vi 

r h_ 
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= J 


Vv 

-Vi 
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-x„ 
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(4) 

( 5 ) 

(6) 


where 


J — x^y v • 


All of the terms involving d ^ and d v are evaluated as finite differences. Equation (6) is the 
matrix form of the metric relations. 
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3. Iterative Interpolation Procedure 


The iterative scheme developed by Benek et al. (ref. 2) is applicable to three- 
dimensional generalized curvilinear coordinates. Trilinear interpolation is used. Yarrow 
and Mehta (report to be published as a NASA Technical Memorandum) utilize tricubic 
interpolants within the same iterative scheme. In this section, the iterative procedure is 
presented for a two-dimensional space. 

Since the grid is uniform and rectangular in computational space, the data can be 
readily interpolated as a function of the generalized coordinates by conventional multivari- 
ate interpolation techniques. The interpolant for the function f can be written as 

f^F(^rj). (7) 

The iterative procedure is used to determine the generalized coordinates (£*,77*) which 
correspond to the physical point (x* ,y*) at which the interpolation is to be performed. 
Since each point in computational space is mapped to only one physical point, we can form 
the following interpolants for the Cartesian coordinates: 

X~X{(,T}) 

y ^ v)- 

For a two-dimensional space, bivariate Newton-Raphson iteration is used to solve for 
(£*,77*) such that 


« 



* 

! 




The value of the function / at the point (x*,y*) is then determined from 


(9) 


r = F(crf). 


( 10 ) 


4. Direct Interpolation Procedure 

The direct interpolation procedure proceeds in three steps. First, required partial 
derivative information in physical space is determined using finite-difference approxima- 
tions to the partial derivatives in computational space. The individual cells are then 
analytically transformed into unit squares by a bilinear transformation. Finally, a bicubic 
spline interpolant is determined within each unit square. 

For a function / known at the grid nodes, we require /*, f y , f X x, f xy , fyy Using finite 
differences in computational space, we can determine ft, f v , ftt, ,fw Applying the 
chain rule to equation (3) yields 
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d X x = Zxxdi "I" Vxxdr, 4" da H - 2 ^xVxd^r, 4* Vx dr,r, 

dxy — Zxyd{ 4- flxydi] 4- ZxZyd{{ 4^ {Zz'Hy 4" r lxZy)d^r, 4" VxVy^VV 

dyy = £yyfy 4" f]yyd v 4* £y d(£ + ^y'fjyd^rj 4 Tj y Orjij • 

Equations (3) and (11) may be assembled as 

dxy = Bd iv , 

where 


B = 


( 11 ) 


Similarly, 

where 


d xy = 
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Therefore B = C 1 . 

In block matrix form, B can be written as 


( 12 ) 


(13) 


B = 


Bi 0 

B 2 b 3 


(14) 


where 



c 



£zx 

i — 
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R* 
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Similarly, 
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(15) 
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where 


Cl = 

x e Vi 

, c 2 = 

x (v 

y«" 

yoi 

, and C3 = 

1 

£ » 

2®*y* 

Vi 2 ' 
ytVv 


. x v y v . 


_ X V V 

Vvv . 


_ x n 2 

V rj 
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Now B — C 1 gives 


B i = G'i 1 (the metric relations), (16) 

B 3 =C 3 -\ (17) 

and B 2 = -C s - 1 C 2 Cr 1 . (18) 

Substituting equations (16) and (17) into equation (18) gives 

B 2 = -B 3 C 2 B 1 . (19) 


Finally, use of the metric relations gives 

£xz = — [?/*; — + 2A;2/£ x ryq ~ x ijVi i Vtt + — a^y^ 2 /t/t? ] 

»?x* = - JJ [-</(y7, 2 z« + 2y ? 2 y, - y* 3 ®,,,, + *<y«, 3 y« - 2 x ( y v y ( y (v + a^y^y,,] 

[ — Vv x v x tz "F yv(vv x t "b — yvVi x i x w + x t] y^ya 

- x v{yv x i + *nVt)yev + ** iVt x tVrm] 

Vxi I = -J3 [ytVv x n x €t - ydVv x i + x vVt) x tv + Vt 2x i x vv ~ HVn^VU 

+ **(y,** + * fl y<)y^ - x^y^w] 

£yy 'j 3 [yv x v x ti ~ 2 y^x^x^x^ 4- y-qX^ x ^ — x ^ y^ + 2x^ x^y^ — x^x^ y^ 

Vyy [~ y( x *i x (( "b 2y{X v X£X{ v V( x ^ x vv ~b x^x^ y^ — 2x r< x^ + a;^ y^j;]. 

(2°) 

Therefore, all of the elements of matrix B can be expressed in terms of the elements 
of matrix C. As a result, d xy f can be determined given d^f. The elements of C and d^,f 
can be approximated using finite differences. An alternative procedure is to calculate f x 
and f y given and from equation (3) and then to use finite differences to determine 
(fx)t,{fy)$Afx)viifv)v The required values of f xx ,f xy , and f yy can then be calculated 
using equation (3) again. The first method involves more computer storage but it also 
involves a smaller stencil. 

Given f x , f y , f xx , f xy , f yy at the corners of a given cell, the problem remains to inter- 
polate / on an arbitrary quadrilateral. We now consider the new two-dimensional space 
given by (p, q) shown in figure 1. The quadrilateral in physical space is related to a unit 
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square in (p, q) space by a bilinear mapping as follows. (Without loss of generality, we will 
assume the lower-left-hand corner to be at (0,0).) 


x = ap + bq + cpq (21a) 

y = dp + eq + fpq (216) 


Note that a different mapping is used for each quadrilateral, or cell, in the grid. This con- 
trasts with the generalized curvilinear coordinate transformation in which a single smooth 
numerically generated mapping is applied to the entire grid. Consequently, analytical 
mappings can be utilized. 



Figure 1. - Two-dimensional space given by (p, q). 

The coefficients of the bilinear mapping are given by 

a = £ 4 , 6 = a? 2 > c = £3 — £4 — £ 2 > 
d = 2 / 4 , e = y 2 , / = J/3 — J/4 — 2/2 • 

The required derivatives in (p, q) space are given (as in equation (12)) by 
Op = XpOx "I" VpOy 

Oq — XqOx "f Vqdy 

dpq = Xpqd x "4* VpqOy d" XpXqO X x d~ (®pj/g d" ^qVp^Oxy d” VpVqOyy i 

where 


( 22 ) 


(23) 


£p = « + cq, y p = d + fq, x q = b + cp, y q = e + /p, x pq = c, y pq = /, 

and we have used x pp = x qq = y pp = y qq = 0. 

Given a location in physical space (xo,t/o)> we must find the corresponding (po,</o)- 
From equation (21a), we have 


9o 


£0 — apo 
6 d~ cp 0 


(24) 


6 



Substituting this into equation (21b) as follows: 


yo = dpo + (e + fpo)( X ,° — ) (25) 

o + cpo 

gives the following quadratic for po 

Po 2 (cd - af) + po{-cy 0 + bd + x 0 f - ae) + (-yob + ex 0 ) = 0. (26) 

The solution of equation (26) gives two values of po- The corresponding values of <70 are 
determined from equation (24). One location (po,</o) will not lie within the unit square 
and hence can be discarded. 

Using fxifyifxxifxytfyy as calculated from equation (12), we can determine / p ,/ g , 
and f pq from equation (23) at each corner of the unit square in (p, q) space. We then 
use the standard formula for bicubic interpolation within a unit square (see Appendix). 
The bicubic is finally evaluated at (po,<7o), as calculated from equations (26) and (24), 
respectively. 


5. Conclusions 

A procedure has been presented for interpolation on a structured, curvilinear, two- 
dimensional grid. In contrast to existing methods (ref. 2), the procedure avoids the need for 
iteration. A piecewise bicubic spline is used, leading to a smooth and accurate interpolant 
which is capable of preserving derivative information generated by a flow solver based on 
finite differences. Potential applications include multiblock algorithms using overlapping 
grids in the solution of problems in computational fluid dynamics, surface definition, and 
graphical representation of data. Extension to three dimensions is straightforward but 
nontrivial. 
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Appendix 

Locally Bicubic Spline Interpolation on a Unit Square 

In this appendix, we give a piecewise bicubic spline formulation for interpolation 
within a unit square. De Boor (ref. 5) presents a technique for determining the required 
partial derivatives at the corners of each cell in a rectangular grid such that C 2 continuity 
is obtained for the domain. In the formulation used here, the required partial derivatives 
at the corners of the unit square are given by centered difference formulas, leading to C 1 
continuity, analogous to the univariate cubic Bessel spline given by de Boor (ref. 6). This 
formulation has the advantage of producing a local interpolant. 

The bicubic spline interpolant on a unit ( p , q) square can be written as 


3 

/(?»?) — ^ ] 7 mnP 9 ? 0 < p, 9 ^ 1? {A ~~ 1) 

m,n=0 


where 

[7m«] = AKA 1 , 
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and 


r/(o,o) /,(o,o) 

TS /p(OjO) fpq( OjO) 

” /(1,0) /,(1,0) 

L/,(i,o) / M (i,o) 


mi) /,( 0,1)] 

/ p ( 0 , 1 ) /„( 0 , 1 ) 

4 ( 1 , 1 ) 
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